%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% APPENDIX: ESTIMATE THE EFFECT OF TRANSPORTATION COSTS ON DEFORESTATION 
% BASED ON A SEMIPARAMETRIC QUANTILE IV MODEL - SPQIV

% "DEFORESTATION IN THE AMAZON:
% A UNIFIED FRAMEWORK FOR ESTIMATION AND POLICY ANALYSIS"

% by Eduardo Souza-Rodrigues

% This version: November 2018

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MODEL:

% G(Y,U) = X'beta(U) - TC

% where:
% Y     = share of deforested area
% TC    = transportation cost to the nearest port
% X     = exogenous regressors
% U     = municipality-level unobservable 
% D     = [D1 D2], Euclidean distances to: 
%         (i) nearest port and (ii) nearest capital
% beta  = finite dimensional parameter
% G     = infinite dimensional parameter such that 0 < G^{-1} < 1
%         (Logit case: G(x) = log(x) - log(1-x) )

% Moment Restrictions:
% E(1[G(Y,u) - X'beta(u) + TC < 0] - u | X, D) = 0

% Estimation method: 
% Penalized Sieve Minimum Distance Estimator (PSMD)

% Unit of observation: municipalitites in the Brazilian Amazon

% OBS: Initial Guess for the minimization routine to estimate SPQIV depends
% on the parametric estimates obtained in "dem_def_reg.m"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (A) SET UP PARAMETERS FOR SIEVES ESTIMATION (PSMD)


% I. Parameters for the Estimation Procedure

% Penalization Terms
lambda1  = 0.000001;  % Penalization for "roughness": || [D2(h)]^2 || < inf
lambda2  = 0;         % Penalization for violations : 0 < h < 1.
lambda3  = 0;         % Penalization for violations of monotonicity of h.

% Naive or smooth version for criterion function:
% approx = 0, use indicator function in the moment restriction
% approx = 1, use smooth (logistic) approximation for indicator function
%             in the moment restriction
approx = 0;

% Coeficient that controls the smoothness of the logistic approximation (if use approx = 1)
p_ap = 1;

% Number of initial values for the optimization routine:
V=5;
%V=1;

% II. Define Approximation Functions

% (a) Approximation function for G(Y): 
% Use Artificial Neural Network, ANN(kh)
% Order
kh = 3;
% Dimension
dimh = 3* kh;
% Total number of parameters to estimate
dim  = dimh + dx;

% (b) Approximation for the Moment Condition Function, m(Z): 
% Use Polinomial Splines, SPLINE(km,qm) (but only for D = [D1 D2])
% Break points
qm = [0.01 0.25 0.75 0.99];
%qm = 0.1:0.1:0.9;
% Order
km = 2;
% Dimension
dimm = km + size(qm,2) - 1;
% Generate P matrix of instruments:
P1   = splibas(quantile(Z(:,1),qm),0,km,Z(:,1));
P2   = splibas(quantile(Z(:,2),qm),0,km,Z(:,2));
P    = [P1 P2 Z(:,1).*Z(:,2) X];
% Dimension of matrix P 
dimp = size(P,2);
    
% Check dimensions
if (dim > dimp)
    error('More parameters than instruments')
end


% III. Initial Guess for Minimization Routine

% (a) Initial guess for beta: normalized IVQR estimates -- Only at Median
beta_guess = -beta_ivqr(2:dx+1,ts(3))./beta_ivqr(1,ts(3)); 
%for i=1:t   % Other quantiles
%    beta_guess(:,i) = -beta_ivqr(2:dx+1,i)./beta_ivqr(1,i);
%end

% (b) Initial guess for G(Y): approximate (normalized) logit function by ANN
% OBS: ANN is characterized by three vectors of parameters: alpha, mu, and sigma

% Evaluating points
y = nodeunif(10001,myzero,1-myzero); 

% Initial Guesses for Minimization Routine (for deterministic approximation)
mu    = [0.25 0.5 0.75]';
sigma = ones(kh,1);

% Normalized Logit Function (to be approximated by ANN
h = (log(y./(1-y)) - beta_ivqr(dxc,ts(3)))./(-beta_ivqr(1,ts(3))); 

% Options for minimization routine
options = optimset('Display','iter');

% Find Best Deterministic Approximation (search over [mu;sigma])
musigma = fminsearch('neural_net',[mu;sigma],options,h,y,kh);
%musigma = fminunc('neural_net',[mu;sigma],options,h,y,kh);

% Store result
mu_guess    = musigma(1:kh,1);
sigma_guess = musigma(kh+1:2*kh,1);

% Store corresponding alphas
[f,alpha_guess,g] = neural_net(musigma,h,y,kh);

% Check deterministic approximation
%plot(y,g*alpha_guess,'o',y,h,'g')

clear h musigma f

% For other quantiles  
%for i=1:t 
%    h(:,i) = (log(y./(1-y)) - beta_ivqr(dxc,i))./(-beta_ivqr(1,i)); 
%    % Options
%    options = optimset('Display','iter');
%    % Find [mu;sigma]
%    musigma(:,i) = fminsearch('neural_net',[mu;sigma],options,h(:,i),y,kh);
%    %musigma(:,i) = fminunc('neural_net',[mu;sigma],options,h(:,i),y,kh);
%    % Store result
%    mu_guess(:,i)    = musigma(1:kh,i);
%    sigma_guess(:,i) = musigma(kh+1:2*kh,i);
%    [f(:,i),alpha_guess(:,i),gi] = neural_net(musigma(:,i),h(:,i),y,kh);
%    g(:,:,i) = gi;
%end
% Check approximation
%if check_ap == 1
%    for i=10:10:90
%        gi = g(:,:,i);
%        i
%        plot(y,gi*alpha_guess(:,i),'o',y,h(:,i),'g')            
%        pause
%    end
%end
%clear h musigma f
    

% (c) Stack parameters with initial guesses
initial_guess = [alpha_guess; mu_guess; sigma_guess; beta_guess];

% Other initial values by perturbing the initial guess
rng('default')  % Set seed
guess = zeros(dim,V);
guess(:,1) = initial_guess;
for v=2:V
     guess(:,v) = initial_guess + randn(dim,1);
end


%% (B) ESTIMATION PROCEDURE

% Define Matrices to store results
alpha_np_v = zeros(kh,V);     % alphas
mu_np_v    = zeros(kh,V);     % mu
sigma_np_v = zeros(kh,V);     % sigma
beta_np_v  = zeros(dx,V);     % betas
q_np_v     = zeros(1,V);      % Criterion Function

% For other quantiles
%alpha_np = zeros(kh,t);     % alphas
%mu_np    = zeros(kh,t);     % mu
%sigma_np = zeros(kh,t);     % sigma
%beta_np  = zeros(dx,t);     % betas
%q_np     = zeros(1,t);      % Criterion Function

% Estimate parameters using PSMD estimator
for v=1:V
    
    % (i) Search over the parameters    
    % Options
    options = optimset('Display','iter');
    
    % Minimization
    [theta,Q,exitflag]= fminsearch('npqiv_ann_log',guess(:,v),options,Yt,TC,X,P,kh,0.5,approx,p_ap);
    
    % Improve results
    while exitflag == 0
        [theta,Q,exitflag]= fminsearch('npqiv_ann_log',theta,options,Yt,TC,X,P,kh,0.5,approx,p_ap);    
    end
    
    % (ii) Store quantile results
    alpha_np_v(:,v) = theta(1:kh,1);
    mu_np_v(:,v)    = theta(kh+1:2*kh,1);
    sigma_np_v(:,v) = theta(2*kh+1:3*kh,1);
    beta_np_v(:,v)  = theta(3*kh+1:dim,1);
    q_np_v(v)       = Q;      % Value of Criterion Function
    
end

% Find Minimum across initial guesses
index = find(min(q_np_v));

% Store results for the minimum across V
alpha_np = alpha_np_v(:,index);
mu_np    = mu_np_v(:,index);
sigma_np = sigma_np_v(:,index);
beta_np  = beta_np_v(:,index);

   
% Fix quantiles - Estimate only at quantiles ts = [10 25 50 75 90] 
%for i=1:size(ts,2)
    
    %disp('Quantile')
    %disp(tss(i))
    
    % (i) Search over the parameters
    
    % Options
    %options = optimset('Display','iter');
    
    % Minimization
    %[theta,Q,exitflag]= fminsearch('npqiv_ann_log',initial_guess(:,ts(i)),options,Yt,TC,X,P,kh,tau(ts(i)),approx,p_ap);
    
    % Improve results
    %while exitflag == 0
        %[theta,Q,exitflag]= fminsearch('npqiv_ann_log',theta,options,Yt,TC,X,P,kh,tau(ts(i)),approx,p_ap);    
    %end
    
    % (ii) Store quantile results
    %alpha_np(:,ts(i)) = theta(1:kh,1);
    %mu_np(:,ts(i))    = theta(kh+1:2*kh,1);
    %sigma_np(:,ts(i)) = theta(2*kh+1:3*kh,1);
    %beta_np(:,ts(i))  = theta(3*kh+1:dim,1);
    %q_np(ts(i))       = Q;      % Value of Criterion Function
    
%end
      

%% (C) SAVE RESULTS OF SPQIV ESTIMATES

if farm == 1
    
    % Estimated Coefficients - Normalized for Logit IVQR
    beta_ivqr_norm_1  = [-1;beta_guess]; 
    beta_spqiv_farm_1 = [-1;beta_np]; 
    
    % Save All
    save beta_ivqr_norm_1   beta_ivqr_norm_1   -ascii
    save beta_spqiv_farm_1  beta_spqiv_farm_1  -ascii
            
elseif farm == 2
    
    % Estimated Coefficients - Normalized for Logit IVQR
    beta_ivqr_norm_2  = [-1;beta_guess]; 
    beta_spqiv_farm_2 = [-1;beta_np]; 
    
    % Save All
    save beta_ivqr_norm_2   beta_ivqr_norm_2   -ascii
    save beta_spqiv_farm_2  beta_spqiv_farm_2  -ascii
    
elseif farm == 3
    
    % Estimated Coefficients - Normalized for Logit IVQR
    beta_ivqr_norm_3  = [-1;beta_guess]; 
    beta_spqiv_farm_3 = [-1;beta_np]; 

    % Save All
    save beta_ivqr_norm_3   beta_ivqr_norm_3   -ascii
    save beta_spqiv_farm_3  beta_spqiv_farm_3  -ascii
        
elseif farm == 4
    
    % Estimated Coefficients - Normalized for Logit IVQR
    beta_ivqr_norm_4  = [-1;beta_guess]; 
    beta_spqiv_farm_4 = [-1;beta_np]; 

    % Save All
    save beta_ivqr_norm_4   beta_ivqr_norm_4   -ascii
    save beta_spqiv_farm_4  beta_spqiv_farm_4  -ascii
    
end

%% (D) SAVE GRAPHS

if results_plot_spqiv == 1
    
    if farm == 2
        
        % Prepare G(Y)
        for j=1:kh
            G(:,j) = normcdf((y-mu_np(j,1)*ones(size(y,1),1))/sigma_np(j,1));
        end
        
        % Save Figure
        spqiv_smallmedium = figure('Visible','off');
        plot(y,G*alpha_np,'b--',y,g*alpha_guess,'g-')
        xlabel('Share of Deforestation','fontsize',16)
        ylabel('G(Y)','fontsize',16)
        title('SMALL-MEDIUM FARMS','fontsize',16)
        legend('SPQIV','Logit IVQR','Location','Southeast')
        saveas(spqiv_smallmedium,'spqiv_smallmedium.png')
        saveas(spqiv_smallmedium,'spqiv_smallmedium.fig','fig')
        
    elseif farm == 3
        
        % Prepare G(Y)
        for j=1:kh
            G(:,j) = normcdf((y-mu_np(j,1)*ones(size(y,1),1))/sigma_np(j,1));
        end
        
        % Save Figure
        spqiv_mediumlarge = figure('Visible','off');
        plot(y,G*alpha_np,'b--',y,g*alpha_guess,'g-')
        xlabel('Share of Deforestation','fontsize',16)
        ylabel('G(Y)','fontsize',16)
        title('MEDIUM-LARGE FARMS','fontsize',16)
        legend('SPQIV','Logit IVQR','Location','Southeast')
        saveas(spqiv_mediumlarge,'spqiv_mediumlarge.png')
        saveas(spqiv_mediumlarge,'spqiv_mediumlarge.fig','fig')
        
    elseif farm == 4
        
        % Prepare G(Y)
        for j=1:kh
            G(:,j) = normcdf((y-mu_np(j,1)*ones(size(y,1),1))/sigma_np(j,1));
        end
        
        % Save Figure
        spqiv_large = figure('Visible','off');
        plot(y,G*alpha_np,'b--',y,g*alpha_guess,'g-')
        xlabel('Share of Deforestation','fontsize',16)
        ylabel('G(Y)','fontsize',16)
        title('LARGE FARMS','fontsize',16)
        legend('SPQIV','Logit IVQR','Location','Southeast')
        saveas(spqiv_large,'spqiv_large.png')
        saveas(spqiv_large,'spqiv_large.fig','fig')
        
    end
    
end
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
